//============================================================================
// Name        : cppmcce.cpp
// Author      : x
// Version     :
// Copyright   : Your copyright notice
// Description : Hello World in C++, Ansi-style
//============================================================================

#include <iostream>
#include <map>
#include <ctime>
#include <cstdio>
#include <cstdlib>
#include "Protein.h"
#include "Conformer.h"
//#include "oppzip.h"


using namespace std;
using namespace mcce;


int main() {
	time_t time_start, time_end;
	time_start = time(NULL);

	// read "run.prm" file
	Protein prot;
	Protein::loadEnv();

	// read "head3.lst" file
	prot.loadConfFromHead3("head3.lst");
	prot.adjustHead3Flag();

	cout << prot.getResidues().size() << " residues in total" << endl;

	// read "energies.opp"
	prot.loadConfOpp("energies.opp", Protein::env.scale_vdw, Protein::env.scale_ele);

	time_end = time(NULL);
	cout << "mc setup time: " << time_end - time_start << endl;

	srand(2);
	for (int iPoint=0; iPoint<Protein::env.titr_steps; iPoint++) {
		float ph = Protein::env.titr_ph0;
		float eh = Protein::env.titr_eh0;
		if (Protein::env.titr_type == "p") ph = Protein::env.titr_ph0 + iPoint * Protein::env.titr_phd;
		else eh = Protein::env.titr_eh0 + iPoint * Protein::env.titr_ehd;


		// annealing
		prot.setupMC();
		prot.setState();

		prot.mcSample(Protein::env.monte_nstart, ph, eh);
		prot.mcSample(Protein::env.monte_neq, ph, eh);

		// reduce conformers
		prot.reduceFreeRes();

		// sampling
		prot.setupMC();
		vector<float> mcRunsOcc(prot.getNTotalConf());
		for (int jrun=0; jrun<Protein::env.monte_runs; jrun++) {
			prot.setState();
			prot.mcSample(Protein::env.monte_nstart, ph, eh);

			prot.mcSample(Protein::env.monte_niter, ph, eh);

			for (int k=0; k<prot.getNTotalConf(); k++) mcRunsOcc[k] += prot.getConflist()[k]->getOcc();
		}
		for (int k=0; k<prot.getNTotalConf(); k++) prot.getConflist()[k]->setOcc(mcRunsOcc[k]/Protein::env.monte_runs);

		prot.saveConfOcc();
	}
	prot.printConfInfo();

	time_end = time(NULL);
	cout << "Total running time: " << (time_end - time_start) << endl;

	return 0;
}
